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Abstract 

Developments in dynamical systems theory provides new support 
for the discretisation of pdes and other microscale systems. By sys- 
tematically resolving subgrid microscale dynamics the new approach 
constructs asymptotically accurate, macroscale closures of discrete 
models of the pde. Here we explore reaction-diffusion problems in 
two spatial dimensions. Centre manifold theory ensures that slow 
manifold, holistic, discretisations exists, are quickly attractive, and 
are systematically approximated. Special coupling of the finite ele- 
ments ensures that the resultant discretisations are consistent with 
the PDE to as high an order as desired. Computer algebra handles 
the enormous algebraic details as seen in the specific application to 
the Ginzburg-Landau equation. However, higher order models in 2D 
appear to require a mixed numerical and algebraic approach that is 
also developed. Being driven by the residuals of the equations, the 
modelling here may be straightforwardly adapted to a wide class of 
reaction-diffusion differential and lattice equations in multiple space 
dimensions. 
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1 Introduction 

Here we extend the dynamical systems 'holistic' approach to the macroscale 
discrete modelling of the class of two dimensional, homogeneous, nonlinear 
reaction-diffusion equations 

^ = V- [f(u,Vu)Vu] +ag(u). (1) 

Our approach systematically models subgrid scale processes with the aim 
of providing an effective closure for the macroscale discretisation. Sections 
[2] and [5] discuss two distinct avenues of theoretical support: that of cen- 
tre manifold theory, and consistency, respectively. Such macroscale closures 
should enable a relatively coarse numerical grid to significantly improve com- 
putational speed and stability. 

As a particular example. Sections [3] and H] explore in some detail the 
real valued, two dimensional, Ginzburg-Landau equation obtained from the 
PDE ([T]) with reaction g = u — and with constant diffusion f = 1 , namely 

^ = V^u + a(u-u^). (2) 
ot 

We choose this 2D real Ginzburg-Landau equation as a prototype reaction- 
diffusion PDE because it is well studied and its dynamics well understood [HI 
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[131 e.g.]. The non-trivial stable and unstable steady states of the 2D Ginzburg- 
Landau equation measure the accuracy and effectiveness of various accu- 
racy models in this application. 

The macroscale discretisation is based upon dividing the domain into 
overlapping finite elements. Following analogues in one dimension [231 HBl 
e.g.], neighbouring elements are coupled with the non-local conditions (jl])-® 
with the strength of the coupling parametrised by y. Such coupling of over- 
lapping elements appear analogous to other multiscale methods [5l[29l[71 e.g.]. 
Section 12.21 then discusses how centre manifold theory assures us of the exis- 
tence of a slow manifold that is an exactly closed discrete model. Further this 
slow manifold discretisation is exponentially quickly attractive. Although we 
cannot find this exact slow manifold, theory asserts it may be approximated 
to any asymptotic order in the strength of the interelement coupling y and 
the nonlinearity cc. The overlapping finite elements together with the special 
coupling conditions (j4j)-(l5|) assure us that the resultant macroscale discrete 
models are also consistent with the dynamics of the reaction-diffusion pde 
(Section ED. 

Section [H] outlines the construction, consistency and predictive accuracy 
of second order asymptotic approximations to the macroscale discretisation 
of the Ginzburg-Landau pde ([2]). To extract another order of accuracy from 
the algebra, we find (for the first time) the adjoint operator of the diffusion 
operator on the elements with the nonlocal coupling conditions. The null 
space of this adjoint, strikingly similar to a Galerkin basis, enables us to 
use an integral solvability condition to construct the third order discrete 
model (fTTj) . 

However, higher order holistic models cannot be found analytically. This 
inability to construct algebraic approximations is one major difference be- 
tween systems in one and multiple spatial dimensions. Section H] explores how 
to numerically construct the subgrid field and its evolution in 2D reaction- 
diffusion equations using the 2D Ginzburg-Landau equation as an example. 
We find that even a relatively coarse subgrid microscale resolution is adequate 
to accurately predict the macroscale dynamics. 

Because of the faithful resolution of subgrid structures and interelement 
interactions, the resulting discrete models are algebraically complicated as 
seen, for example, in the third order model ( fT7|) . Thus users may prefer, as 
they often do now, to use models of the nonlinear dynamics of lower order. 
Then higher order discretisations derived via this approach provide good local 
estimates of the local error in a lower order simulation as it is computed on 
the fly. 

A further application should be to the 'equation-free patch' methodology 
for simulating multiscale systems [HI [291 e.g.]. This work suggests, analogous 
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Figure 1: The discretisation of a 2D domain into square elements: The 
t, jth element (solid) is centred upon the grid point {xi,yj)] Eij overlaps 
neighbouring domains to extend to the neighbouring grid points. 

to one dimensional dynamics [2^ [53] , that straightforward coupling condi- 
tions will also empower efficient and accurate patch dynamics in multiple 
space dimensions. 

2 Divide the domain into square elements 

We place the discrete modelling of two dimensional, reaction-diffusion equa- 
tions within the purview of centre manifold theory by dividing the domain 
into overlapping square elements, as shown schematically in Figure [H The 
discretisation is similar to that for shear dispersion in a long thin channel [T3] . 
The significant difference here is that the discretisation of elements is in both 
spatial dimensions, not just along the channel as for the shear dispersion ap- 
plication. 

2.1 Extend the non-local IBCs to 2D 

In this initial study, for two dimensional, reaction-diffusion equations we 
divide the domain into a set of overlapping square elements. Figure [H Define 
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a grid of points (xt,ijj) with, for simplicity, constant spacing h. The t, jth 
element, Etj? is then centred upon [Xi^y^] and of width Ax = Ay = 2h. Let 
Vij(x, ijit) denote the field in the i, jth element and so evolves according to 
the reaction-diffusion pde ([T]); that is, 

^ - V ■ [fKj, Vv,j)Vvtj] + (xgKj) . (3) 
The original field u[x,y,t) is then prescribed by Vij[x,y,t] when [x,y] G 

The evolution of the field over the whole domain then depends upon how 
the elements are coupled together. To couple the dynamics of each element to 
its neighbours we use coupling 'internal boundary conditions' (iBCs) around 
the t, jth element of 

Vi,i(xi±i,-y,t) =TVi±ij(xi±i,-y,t) + (1 -y)vij(xi,-y,t), \y -y^\ <h, 

(4) 

Vi,][x,yj±^,t) =yvij±i(x,'yj±i,t) + (1 - yjvy (x.ijj, t), |x - Xi| < h. 

(5) 

These iBCs are a natural extension to 2D of iBCs established for ID dynam- 
ics [231 e.g.]. The crucial feature is: with y = the elements are effectively 
isolated from each other, dividing the domain into decoupled elements with 
consequently independent dynamics; whereas with y = 1 these iBCs ensure 
sufficient continuity between elements to recover the original problem over all 
space. Such coupling of overlapping elements appears analogous to the 'bor- 
der regions' of the heterogeneous multiscale method [5], e.g.], to the 'buffers' 
of the gap-tooth scheme [29l e.g.], and to the overlapping domain decomposi- 
tion that improves convergence in waveform relaxation of parabolic pdes [71 
e.g.]. 

For definiteness in theoretical support, let there be m elements in the 
domain with the field required to be periodic in both x and y. For example, 
the elements may form a ^/^r^ x ^/^r^ grid in the domain (any factorisation 
of m is feasible). In principle, the 2D elements could be any shapes, regular 
or irregular; square elements appear to be easiest to start with. In this initial 
work we also avoid issues of physical domain boundary conditions on the 2D 
domain by adopting doubly periodic solutions. (Physical domain boundary 
conditions have been explored for ID domains [271 CSl 6.g.].) These adoptions 
enable straightforward theoretical statements of support. 
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2.2 Centre manifold theory supports discrete models 

This section describes in detail how the iBCs (l4])-(l5l) lead to centre manifold 
theory supporting an accurately closed, discrete model for react ion- diffusion 
systems ([I]) via its dynamics ([3]) on overlapping elements. 

A homotopy in the coupling parameter y connects the physically relevant 
discretisation to a tractable base. When parameters a = y = both the 
reaction and the coupling on the right-hand side of the iBCs (jl])-© disappear. 
The elements are then effectively isolated from each other and so the resultant 
diffusion in the PDE ([3]) is particularly simple: exponentially quickly in time, 
the field Vt j becomes independently constant within each element. We use 
this family of piecewise constant solutions as a basis for analysing the case 
when the elements are coupled together, y ^ . Particularly interesting is 
the approximation for full coupling, y = 1 , when the PDE ([T]) is effectively 
restored over the whole domain because iBCs (HJ-Q then ensure sufficient 
continuity between adjacent elements as described previously for ID pdes [28l 
[IS1I2SI e.g.]. 

The support of centre manifold theory is based upon a linear picture of 
the dynamics. Adjoin the dynamically trivial equations 

and consider the reaction-diffusion dynamics in the extended state space 
(vij(x,'y],y, a). In this extended state space, points oc — 'y — and Vij = 
constant are equilibria of the diffusion hence these form a subspace of 
equilibria, Eq — {(vtj,0,0)}, in the extended state space. Linearized about 
each of the equilibria in Eq, the PDE for perturbations v[^[x,y,t) within each 
element is then the constant coefficient diffusion PDE 

= fijV for {x,y] e U,j for each i, j, (7) 

where the constant diffusivities fij = f(vtj,0). These pdes are decoupled 
because they are to be solved with the y = iBCs 

v(j(x,-yj±i,t) = v(^j(x,-yj,t], |x-Xi|<H. (8) 
Thus the following linear eigenmodes are associated with each element: 
cx = y = , 

v{j oc e^^'i''^''* X cos [k7r(x - Xi_i/2)/H] x cos [ln[y - Yj_i/2]/H] , 
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inside the t, jth element for all integers k, I = 0, 1 , 2, . . ., where the decay rate 
of each mode is 

Aij,k,i = -fi,j 

together with the two trivial modes that firstly y = constant and a — — 0, 
and that secondly ex. — constant and y — v[- = . In a spatial domain 
with m elements and when all diffusivities fij > , then all eigenvalues are 
negative, — ftjTr^/h,^ or less, except for m + 2 zero eigenvalues. Of the m + 2 
zero eigenvalues, one is associated with each of the m elements and two come 
from the trivial equations for the parameters. That is, the slow subspace 
is {(vij,y, (x)} for constant V|j. The above arguments establish the following 
corollary of a centre manifold existence theorem ([3l p. 281] or [31], p.96]). 

Corollary 1 (Existence) Provided the nonlinear diffusivity f and reaction g 
in ([3j) are sufficiently smooth, and all fij > then a m + 2 dimensional slow 
manifold M. exists for ([3]) -([6]) in some finite neighbourhood of the subspace Eq 
of equilibria^ 

The slow manifold M. is parametrized both by the two parameters y and a, 
and by a measure of the field in each element; we use the grid value Utj(t) = 
Vij^XijUj, t) to measure the field in the i, jth element. Using u to denote the 
vector of such parameters, we write the slow manifold M. as 

Vij =vtj(x,'y;u,y,a). (10) 

These functions Vt j(x,'y; u,y, a), that Sections [3] and H] construct for the 
Ginzburg-Landau pde ([2]), resolve the subgrid scale physical structures as a 
function of the neighbouring grid values in u. On this slow manifold M. the 
grid values Utj evolve deterministically 

duij/dt = uy = gij(u,y, a) , (11) 

where gtj- is the restriction of (l3])-(l6l) to the slow manifold M.. In essence, 
this closure of the grid scale dynamics comes from the accurate resolution of 
the subgrid scale structures. 

Using the value of the field at the grid points to parametrise the slow 
manifold provides the necessary 'amplitude conditions' to close the problem: 

Uij =v(xi,-yj;u,y,(x). (12) 

^Keep clear the distinction between centre manifold theory and the slow manifolds 
discussed here: the theory applies to systems where the real part of the eigenvalues of 
critical modes are zero; whereas here wc explore and construct slow manifolds because 
here the eigenvalues arc precisely zero. 
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Many other amplitude conditions are possible such as defining the 'ampli- 
tudes' Uij to be the mean field over the i, jth element. However, using the 
grid values are simple, traditional, and have a direct physical meaning. 

Centre manifold theorems [3], [SU e.g.] also support the following crucial 
emergence and approximation properties. 

Corollary 2 (Emergence and approximation) Provided the nonlinear dif- 
fusivity f and reaction g in are sufficiently smooth, and all ftj > , then 

• every solution of the reaction- diffusion dynamics ©-([H]) that stays 
within a neighbourhood of the slow manifold Ai, ( ITOi) . approaches ex- 
ponentially quickly a solution of the discrete model ( ITT]) on the slow 
manifold (fTOj) .- and 

• the order of error in asymptotically approximating the slow manifold 
and its evolution, ( fTOl) -(fTTl) . is the same as the order of residuals of 
the governing equations ([3]) -([6]). In particular, because the base equi- 
libria form a subspace Kq, here the approximation is global in the grid 
values Uij, it is local only in the two parameters y and a. 

3 A slow manifold discretisation of the Ginzburg— 
Landau equation 

We now explore a slow manifold discrete model for a specific reaction-diffusion 
system, the 2D Ginzburg-Landau equation equation ([2]). Substituting ( fTOl) 
and (ITT!) , the PDE we solve to form the model is obtained by equating the 
PDE for dvi j/dt to that obtained by the chain rule: 

^-L|5i9M-VV, + «K,-vf,,). (13) 

To construct the slow manifold f|TOl) - f|TTl) by solving the PDE f|T3|) with cou- 
pling and amplitude conditions involves considerable algebraic detail. These 
algebraic details of the construction are handled straightforwardly by itera- 
tion in computer algebra [261 e.g.]. The specific procedure used here, docu- 
mented elsewhere [18], solves the equations using iteration to drive to zero 
the residuals of the governing differential equation (fT3ll and its interelement 
coupling iBCs (II])-©. Hence by the Approximation Corollary [21 we construct 
correspondingly accurate approximations to the slow manifold of (1131) . These 
approximations, upon setting coupling parameter y = 1 , form 2D discrete 
models of the the 2D Ginzburg-Landau PDE ([2]). 
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One consequence of using computer algebra is that there is no need to 
record in this article most of the considerable algebraic detail in constructing 
the models. Those wishing to verify the correctness of the results recorded 
herein should download and examine the corresponding technical report [18] 
that details the precise computer algebra procedure. Because the algorithm is 
based upon driving the residuals to zero, the critical aspect of the procedure is 
simply the correct coding of the computation of the residuals of the governing 
equations. 

The 0(y^ + oc') holistic discretisation Satisfying the pde and iBCs 
to residuals of 0(y^ + oc') the computer algebra procedure [IHl §2.2] gives 
subgrid fields which are too complex to record here. The corresponding 
evolution of the grid values on the slow manifold are 

where the centred difference operator applies in both spatial dimensions, 
6 = Ui+i J + U|_i J + Uij+i + Uij_i - 4uij , 

- 4(ui+i J + Ui_i J + Uij+i + Uij_i ] + 1 2uij , 

The model f ll4p is simply the extension to two spatial dimensions of the 
0(y^ + <^^) holistic model of the ID Ginzburg-Landau equation [T7j . 

The holistic discrete model has the dual justification of consistency with 
the PDE in addition to the justification provided by centre manifold theory. 
As proven in Section [5l consistency for such discrete models follows from 
the coupling iBCs (H])-© [23]. Set the coupling parameter y = 1 in the 
discrete equation ffT^ to recover the holistic discrete model of the Ginzburg- 
Landau PDE (12]) in 2D. To test consistency, we expand the finite differences 
of (|T4l) in a Taylor series in the grid spacing h [181 §2-5] to find the equivalent 
continuum PDE for the O (y^ + (X^) holistic model ( |T4l) is 

- = vWa(u-u^) + — ulVul^-- ( 3^ + 3^ j +0(«^+h'^) . (15) 

The 0{y^ + <xF) holistic model is 0{\\^ + oc^) consistent, maintaining in 2D 
the dual justification of holistic discretisation found for ID pdes [281 123] . 
Section 15] proves this consistency in some generality for 2D pdes. 
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Figure 2: an example of the subgrid field for the 0(y^, a^) holistic model 
with 4x4 elements on [0, tt] x [0, n] at nonlinearity a = 6 

3.1 Illustration of the subgrid field in 2D 

Here we plot the subgrid fields for a coarse grid solution of the O (y^, a^) holis- 
tic model (obtained from IHM by omitting the term and then evaluating 
at y = 1) of the 2D Ginzburg-Landau equation ([2]). We restrict attention to 
a doubly odd symmetric solution that is 27r-doubly periodic. That is, 

u[x,y,t) =u(x + 27r,-y,t), u[x,y,t) = -u(27r - x,-y , t) 
u(x,-y,t) = ulx.-y +27r,t), u(x,-y, t) = -u(x, 27r - ij, t) . (16) 

Figure [2] shows the subgrid fields for the 0(y^, a^) holistic model with 4x4 
elements on [0,tc] x [0,7r] at nonlinearity cx = 6. The subgrid fields exhibit 
the nonlinear subgrid structure of the 2D Ginzburg-Landau equation and 
its interaction through the iBCs. The subgrid fields are comprised of actual 
solutions, albeit approximate, of the 2D Ginzburg-Landau pde. 

Note the subgrid fields have noticeable jumps across the boundaries of 
the elements. Higher order holistic models should reduce these jumps across 
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Accurate bifurcation diagram for the 2D RGL 
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Figure 3: Accurate bifurcation diagram of the 2D Ginzburg-Landau equation 
system for < a, < 30. It is constructed with a fourth order centered 
difference approximation with 24 x 24 points on [0,n] x [0,n]. Blue curves 
indicate stable steady state solutions and red curves indicate unstable steady 
state solutions. The open squares indicate steady state bifurcations. 

the boundaries as seen for the holistic models of the Kuramoto-Sivashinsky 
equation [TB] . 

The C(y2, a^) holist ic model needs improving We investigate the per- 
formance of the 0(y^, a^) holistic model on coarse grids by comparing its 
bifurcation diagram to an accurate solution. Again we restrict our attention 
to doubly odd symmetric solutions that are 27T-doubly periodic ( fT6|) . 

The bifurcation information is calculated using the continuation software 
AUTO [1] and XPPAUT [U] as outlined for the Kuramoto-Sivashinsky equation 
in MacKenzie's PhD dissertation [TTj. In such bifurcation diagrams the blue 
curves indicate stable steady state solutions and red curves indicate unstable 
steady state solutions. The open squares indicate steady state bifurcations. 

Figure [3] shows an accurate bifurcation diagram of the 2D Ginzburg- 
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Centered 2nd order, 8x8 points 
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Figure 4: Bifurcation diagrams for the 2D Ginzburg-Landau equation with 
8x8 elements on [0,n] x [0,n] for (a) (9(y^, a^) hohstic model (b) second 
order centered difference. The accurate bifurcation diagram is shown in grey. 

Landau pde. It is constructed with a fourth order centered difference approx- 
imation with 24 X 24 points on [0,n] x [0,7r]. The trivial solution undergoes 
steady state bifurcations at cc = 2,8, 18 leading to the unimodal, bimodal 
and trimodal branches respectively. For 1 < a < 30, only the unimodal 
branch is stable and all other branches are unstable. 

Figure m shows a comparison of the 0(y^, oc^) holistic model and a sec- 
ond order explicit centered difference approximation with 8x8 elements 
on [0,7r] X [0,n]. The accurate bifurcation diagram is shown in grey. The 
(9(y^, oc^) holistic model does not perform as well as the second order cen- 
tered difference approximation; this is similar to performance observed for 
the one dimensional Ginzburg-Landau pde [TTj. 

Higher order models need numerical construction To improve the 
accuracy of the holistic discretisation we need to compute higher orders in 
either coupling y, or nonlinearity cc, or both. Improved accuracy occurs 
at higher order in comparable ID problems [16]. However, apparently it 
is not possible to analytically construct higher order subgrid fields in 2D: 
apparently the subgrid fields required for our closures are no longer in the 
class of multivariate polynomials. Instead, numerical methods must be used 
to find the subgrid fields as described in Section HI However, the well known 
'solvability condition' in asymptotic mathematical methods empowers us to 
derive the next order in the evolution, analytically from the residuals, without 
needing to find the next order of the subgrid fields. 
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Figure 5: each element is effectively divided into four subregions by the non- 
locality of the boundary conditions ([H]). To derive the adjoint, label the edges 
of these four subregions as shown. 



3.2 The adjoint provides an extra order of accuracy 



We scrounge an extra order of accuracy from the 'solvability condition' [211 
e.g.] applied to residuals of the next asymptotic order. Because the linear 
operator used to find corrections to the subgrid field is singular — the operator 
necessarily has homogeneous solutions that compose the slow subspace Eq — 
the Fredholm alternative is that the 'right-hand side' of the equation for the 
subgrid fields must lie in the range of the singular operator. This solvability 
condition is enough to determine an extra correction to the evolution. 

Recall from linear algebra that to be in the range of the operator, the 
solvability condition is that the right-hand side must be orthogonal to the 
null space of the adjoint operator. Thus the first task of this section is to 
find the adjoint operator of the linear constant diffusion pde ([71) with its 
boundary conditions ([H]). Second, we find a basis for the null space. Lastly, 
more computer algebra provides the required extra order in the evolution. 

The decoupling of the elements, provided by y = in the boundary 
conditions (IHl), simplifies finding the adjoint: we need only consider each 
element in isolation. Thus define the inner product to be the integral over 
the 1, jth element: 



(v,w) 



vw dx dy 



i + l 



Xt+1 



VW dx dy 



j-1 ■J 



To find the adjoint recognise that each element is subdivided into four sub- 
regions, shown schematically in Figure by the coupling of the boundary 
values to internal values by the boundary conditions ([H]). In addition, there 
exists previously implicit conditions that the subgrid field v and its gradient 
are continuous throughout the element. Then integration by parts, or the 



13 



divergence theorem, transforms the inner product 



V^v,w) = (v, V^w) 



v^w — vwx dy 



+ 



v^yv — vWx dy 



— vWx dy 



M4 



+ analogous integrals on B, C and T, 



where specific parts of the boundary integrals are labelled as shown in Fig- 
ure [5l Using superscripts to denote evaluation, continuity requires v'^''^ — 
and v^"*^ = , and the boundary conditions (JHl) imply v'- = v*^ = v 
the inner product 



^--^ Thus 



(V^v,w) = (v, V^w) + 



'i-l 



dy 



+ (analogous x integrals of y derivatives) 



+ v^(w^-w^-wf- + wf+) dy 

+ (analogous x integrals of y derivatives). 

For the adjoint, these integrals on the right-hand side must vanish for all 
smooth fields v. Consequently, the null space of the adjoint operator satisfies 
Laplace's equation V^w = with conditions: firstly, that w is zero around 
the edges L, R, B and T of the element; secondly, that w is continuous on the 
interior partitions M and C (but its gradients may be discontinuous there); 
thirdly, that — + — w J = ; and lastly, that —w^^ + Wy^ — 

Because of these conditions, the null space of the adjoint is spanned by 
the 'pyramid' w = (1 — |x — Xt|/h,)(1 — |y — Yi|/h.) as displayed in Figure El 
The solvability condition is then that the integral of the subgrid residuals of 
the PDE with this w over the i, jth element determines a correction to the 
model evolution. 

It will not escape your notice that the solvability condition integral paral- 
lels integrals in the Galerkin finite element method. Thus the Galerkin finite 
element method may be viewed as a leading approximation to our systematic 
slow manifold closure of discrete modelling. 
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Figure 6: basis 'pyramid' for the null space of the adjoint operator on an 
element: w = (1 - |£,|)(1 - |ti|) = (1 - |x - Xi|/h] (1 - |-y - Yt|/h]. 

Return to the discrete Ginzburg— Landau model Computer algebra 
readily computes the subgrid residuals of the Ginzburg-Landau pde ((21) to 
the next higher order (detailed elsewhere [IHl §2.3]). Taking the inner product 
with the adjoint null vector w, and remembering contributions from the inter- 
element coupling conditions (|ll)-(l5]), gives the discrete model, in gory detail, 

2 3 

2 

+ }^ [u?,^(2225Vj + 245 Vj - 36^5 Jutj) 

+ utj(-1025Vj + 36{6Vj}" + 6{6^Utj]{6yUt,j} - 144{^i5utj}") 

- 6{^y6y6^Uij}{|j.y5yuJj} - 6{|a.xSxSyUij}{M,x6xuJj} 

+ 12{^i6uJ^}{^5u,j)+ 12{^t6V,j}{^i5uy- |{6Mj«S^Sj^tj} 

+ 3{5 V,j}{5 Vj} - 3{6^uJ^]{6jutj} - 3{62uf,j}{6^ut,j} 

+ 9{5 VjWS Vj} - 85Xj - 66 Vj + ^l^Ki 

+ 0(y4 + a4), (17) 
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This complicated macroscale discrete closure arises from resolving subgrid 
scale nonlinear dynamics within the finite elements. 

To independently check the above model, consider the small element limit 
as h, ^ . Upon setting the coupling parameter y = 1 we expect to recover 
classic consistency. Straightforward computer algebra [181 §2-5] finds that for 
small grid size h. the equivalent partial differential equation to the discrete 
model (fT7|) is 



As found for ID problems with analogous element coupling [23], this slow 
manifold discrete model is consistent to the 2D Ginzburg-Landau equa- 
tion ([2]) to high order in the element size h, both for the linear diffusion 
and the nonlinear reaction. Section [5] proves our coupling conditions con- 
struct consistent discrete models in general. But remember that Corollary [2] 
independently provides strong support for the relevance of the model ( fT71) at 
the finite element sizes used in simulations. 

However, the discrete model fll7p is as high an order of accuracy as we 
can construct analytically. The next Section H] shows how to numerically 
solve for the subgrid scale field in order to construct the macroscale discrete 
model. It serves as a proof of principle for applying the holistic method to 
PDEs of two or more spatial dimensions. 

4 Generally compute 2D subgrid fields nu- 
merically 

The holistic discretisation of PDEs is based upon centre manifold theory [21 
El [121 e.g.] to resolve the subgrid fields and hence more accurately close the 
macroscale discrete model. Here the subgrid field is constructed numerically 
for the Ginzburg-Landau equation ([2]) in 2D. 

New complexities arise. Although the spatial structure is obtained nu- 
merically, the slow manifold, subgrid field is also parametrised by the grid 
values Utj, the interelement coupling parameter y, and the nonlinear pa- 
rameter a. Therefore, the construction involves symbolic parameters. The 
algorithm required to develop the holistic model must efficiently solve the 
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corresponding mixed numerical and symbolic problem. The focus of this sec- 
tion is on this novel numerical construction of the subgrid fields and not the 
performance of the holistic models. 

Numerical construction of the subgrid field introduces errors which are 
separate from the orders of errors of the holistic model. These errors from the 
numerical construction of the holistic discretisation are the major concern. 
The numerical construction of the subgrid field and its evolution has chal- 
lenging details: §4.11 How should the subgrid problem be solved? §4.21 What 
subgrid resolutions will accurately reproduce the analytical holistic models? 
§4.31 What is an efficient implementation? 



4.1 Outline the numerical slow manifold in 2D 

In each element we discretise the microscale subgrid as shown in Figure [71 
At each subgrid grid point we seek the evolution of Uij^k,£ for subgrid indices 
|k|, |£| < n where, for example, Figure [7] shows n = 4. The subgrid is shown 
solid (blue) for this particular example of a 4 x 4 interval subgrid. The 
subgrid field extends to Xi±i and Yj±-[ in order to allow the application of the 
2D non-local iBCs (Ill)-(l5]). The subgrid field does not extend to the extreme 
corners because the subgrid discrete Laplacian applied to the interior points 
does not involve the subgrid field at the corners. 

After discretising the subgrid field in the 1, jth element, classic finite 
differences approximate the spatial derivatives of the subgrid field of the 
Ginzburg-Landau PDE ([2]): 



+ oc(vij,k,«-vf_j_^_«). (19) 

These microscale discretised equations are solved at each of the subgrid points 
inside each of the elements. The elements are coupled with iBCs analogous 
to namely 

Vv j,±n,« = yVi±ij^o,(! + ( 1 - y )ViJ,0,(; 1^1 < TL , (20) 

Vij,ic,±n = yVt,j±i,k,o+ (1 -y)vtj,k,o |k| < n. (21) 

The same centre manifold theorems apply to the system (|T9l)-(l2T]) to assure 
us of the existence, relevance and construction of a slow manifold, macroscale 
discrete model of the dynamics. The macroscale slow manifold to construct 
is that the subgrid field Vij_k,£ — Vij,k,«(u-,y, oc] where, defining uij = Vij^o.o , 
the macroscale grid values u evolve according to Uij = Qi,]['^,y, oc). 
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Figure 7: Example of the 4x4 interval subgrid in 2D; note that such a 
subgrid labelled as "4 x 4" actually extends to be 9 x 9 when overlapped 
with neighbouring elements as shown. 
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We employ an iteration scheme to find the microscale subgrid field and the 
macroscale slow evolution. The initial approximation is that of a constant 
field in each element: Vt j i^^^ ~ Uij such that Uij = gij ~ . Given any 
current approximation, Vij and gtj, we seek an improved approximation 

:= Vij + v(j and := gij + g(j where v(j and g(j are corrections to 
be found in each iteration. At each iteration, the following linear equations 
driven by the current residuals are solved for the corrections v( j and g( ^ 

v{,^,o,o = 0. (22) 

The iteration repeats until all residuals are zero to a specified order of error. 
This iteration scheme follows that for the analytic construction of the subgrid 
field and is documented in full detail in a separate technical report [IHl §3]. 
Centre manifold theory then assures us that the resultant macroscale discrete 
model Uij = gi j(u,y, cx) is accurate to the same order of error. 

For example, for the coarsest possible subgrid field, with just n = 2 
subintervals, the numerical discretisation of the subgrid field gives a low 
order, macroscale model as 

" tI^^'"'^'' + '^^ ^^^^'^''^ ~ +0{y' + oc') , (23) 

Compare with the analytic macroscale discrete model ffT^ : the terms in the 
first line are identical; the same higher order terms appear in the second line 
but the coefficients are in error by 25%. This correspondence is promising 
for such a coarse microscale subgrid discretisation. 

4.2 Low resolution subgrids are accurate in 2D 

How do macroscale models constructed via a numerical microscale, such 
as ( l23l) . compare with analytic macroscale models? We compare in two ways: 
one via the convergence of the coefficients; and the other by the accuracy of 
the predicted bifurcation diagrams. It appears that the microscale subgrid 
need not be of high resolution. 

First look at the coefficients of the (9(y^+a^) models such as (l23l) . Recall 
that the number of microscale subgrid points, from one macro-grid point to 
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Table 1: coefficients in the 0(y^ + a^) models, such as fl23|] . evidently con- 
verge to the correct analytic coefficients, in ffT^ and labelled oo in the table, 
with errors (9(l/rL^) as the resolution of the microscale grid improves. 



n 






ayu?^6^Utj 
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1 

16 


1 

16 


3 

"16 
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5 

64 


5 

64 


15 
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21 

256 


21 

256 


63 
256 


oo 


1 

12 


1 

12 


1 

4 



Table 2: maximum errors in the coefficients of the 0{y^ + (xf) model ( |T71) 
when approximated numerically at three different subgrid resolutions. The 
decrease by at least a factor of four, upon doubling n, indicates quadratic 
convergence. 

n y^/h^ ycx y^/h^ y^a yoc^h^ 
2 0.021 0:062 0:0033 OTR 0:00T6 
4 0.0052 0.016 0.00086 0.040 0.000098 
8 0.0013 0.0039 0.00022 0.010 0.0000061 



the next, in each dimension, is n. The coefficients linear in the coupling 
parameter y and nonlinearity a are exact for all n > 2 , only the higher 
order coefficients vary with subgrid resolution. Table [T] tabulates coefficients 
in these nonlinear terms, those of 0(y^ + oc^) in models such as ( |23|) . for 
some values of n. Evidently the coefficients converge to the exact values with 
error (9(l/rL^). We expect such quadratic convergence from the quadratic 
modelling in fll9l) of the subgrid scale dynamics. 

Similarly, we compare numerically obtained coefficients for the O{y^+0(.^^ 
model (fTTll . Because of the complexity of the model we make a limited 
comparison: for each order in a and y in ( |T7j) . Table [2] reports the largest 
error in the numerically obtained coefficient for three different subgrid scale 
resolutions. Evidently, these maximum errors decrease like to conffim 
the accuracy of the numerical description of the subgrid scale dynamics. 

Incidentally, this agreement between the numerical model, obtained by 
resolving the subgrid scale dynamics at order O (y^ + (X^) , and the analytic 
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(b) 2D-RGL 0(7^, a^), 4x4 interval subgrid 
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(d) 2D-RGL 0(Y^,a^), 8x8 interval subgrid 
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Figure 8: Bifurcation diagrams of the 0(y^, a^) holistic models of the 2D 
Ginzburg-Landau system with 8x8 macroscale elements on [0,n] x [0, n] for 
subgrid resolutions of (a) 2 x 2, (b) 4 x 4, (c) 6 x 6 and (d) 8 x 8 intervals. 
The bifurcation diagram for the analytically constructed model is shown in 
green. 



model (I17p supports the derivation in Section lX^ of an extra order of accuracy 
in the macroscale model without necessarily having to resolve the subgrid 
scale dynamics at the same order. 

Second, we turn to the bifurcation diagram to see the sort of errors in- 
curred in using the approximate models. Figure M shows the bifurcation 
diagrams for the C?(y^, oc^) holistic model of the Ginzburg-Landau system 
for four subgrid resolutions. Here the equilibria shown in green are not the 
accurate solution of the Ginzburg-Landau system, but rather the equilib- 
ria of the analytic 0(y^,<x^) holistic model in 2D (obtained from ( fT4l) by 
omitting the term). Observe that with a subgrid resolution of just 4x4 
intervals the bifurcation diagram for the numerically constructed 0{y^, cc^) 
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2D-RGL 0(/,a^) holistic, 8x8 elements 2D-RGL 0(7^, a^) holistic, 8x8 elements 




Figure 9: Bifurcation diagrams of the (a) C(y^, a^) and (b) C(y^,cx^) 
holistic models of the 2D Ginzburg-Landau system with 8x8 elements on 
[0, tt] X [0, 7r]. The accurate bifurcation diagram without stability information 
is shown in grey. 

holistic model is almost indiscernible from the analytic model over nonlinear- 
ity < a < 20 . Higher subgrid resolutions are indiscernible to even larger 
nonlinearity a. 

As a last comparison of bifurcation diagrams, Figure [9] shows one example 
confirming that by computing to higher order in the interelement interac- 
tions, and resolving the subgrid scale structures numerically, the predictions 
of the numerically derived models do improve. 

Numerically resolving the microscale subgrid structures does generate 
usefully accurate, slow manifold, macroscale discretisations. 

4.3 An efficient computer algebra approach is crucial 

The difficulty associated with the numerical construction of the subgrid field 
is the mixed discrete and symbolic nature of the equations involved in the 
iteration scheme. The size of the system of equations increases as the subgrid 
resolution improves, and the complexity of the symbolic nonlinear residuals 
increases quickly as higher order holistic models are constructed. 

Computer algebra packages such as reduce [2j or mathematica [32] 
have general routines that solve systems of equations such as (1221) . However, 
these solve routines were inefficient for the many equations and complicated 
expressions involved with better subgrid resolution and higher order holistic 
models. Even a low resolution, n — 2 interval subgrid, took many minutes 
in both REDUCE and mathematica using their built in solve routines. 
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Table 3: reduce and MATHEMATICA computational times for numerical 
construction of 0(y^, a^) holistic models of the one dimensional Ginzburg- 
Landau equation for various subgrid scale resolutions, n. 



TL 


REDUCE 


MATHEMATICA 


2 


1.1 S 


70.2 s 


4 


3.1 S 


215.4s 


8 


8.3 s 


367.6 s 


16 


23.7s 


517.7s 



Instead we develop an approach that is practical for implementation with a 
large number of complicated symbolic terms. 

Transform to constant coefficient Recall that at each step of the iter- 
ation scheme we solve a problem for updates to the subgrid field and its 
evolution g( •. Multiply the first (field) equation in fl22|) by h,^ and replace g(j 
by Q' = H^Qt j • Then the left-hand side of the new form of the equations 
has numerical constant coefficients; algebraic expressions only occur in the 
right-hand side. 

Further, the left-hand side of the new equations remain the same for ev- 
ery iteration. Consequently, the first iteration constructs an LU factorisation 
of the left-hand side, which is then used to solve equation ( 1221) for updates 
in every iteration [181 §3]- The LU decomposition is performed once and 
requires approximately ^N^ operations [221 e.g.]. Here the number of equa- 
tions for the subgrid structure are N — (2n + 1 )^ -|- 1 ; for example, N = 290 
for TL = 8 microscale intervals in the subgrid. At each step of the iteration 
scheme the LU factorisation algorithm operates on the symbolic residual vec- 
tor. Perhaps 2D and 3D problems could be solved more efficiently through 
iterative multigrid [T^[T] or incomplete LU factorisation and Krylov subspace 
methods [TUl ED] • Such alternatives remain for later exploration. 

Reduce was faster Computational comparison experiments found that 
the computer algebra reduce was an order of magnitude faster than MATH- 
EMATICA. Table [3] lists the computational time for the reduce and the 
MATHEMATICA implementation for constructing O(y'^,oc^^ holistic models 
of the one dimensional Ginzburg-Landau equation with subgrid resolutions 
of 2, 4, 8 and 16 subgrid intervals. These times were observed on a Pen- 
tium III, 750 MHz processor, with 256Mb RAM, running reduce 3.7, under 
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Windows XP. Table [3] shows the reduce implementation was 20-70 times 
faster than the mathematica implementation (despite the repeated help of 
the mathematica news group). Thus we used Reduce. 

REDUCE is an interpreted language and as such it interprets each com- 
mand and allocate memory as each command is executed at runtime. It 
would be possible to write a purpose built compiled program in some main- 
stream language to study any specific pde. This option is not considered in 
the scope of this article. 

5 Non-local coupling conditions enforce con- 
sistency 

Recall that the constructed holistic models of the Ginzburg-Landau dy- 
namics are consistent with the PDE as the grid size h — > , see equations 
fllSp and fllSp in Section [31 Now we prove that general consistency follows 
from the specific choice of nonlocal interelement coupling conditions (jl])-®. 

We start with a similar theorem to one previously proved for the con- 
sistency of holistic discretisation in one space dimension [23]. The critical 
difference here is in the proof: previously the proof was constructive whereas 
here it is not. Avoiding a constructive proof has two consequences: it is es- 
sential here as we do not know analytic forms for the slow manifold subgrid 
fields in 2D; and the new proof easily caters for nonlinear reaction. Because 
the theorem here is more powerful, an immediate corollary proves consistency 
of a 2D holistic discretisation with the 2D linear PDE. 

Theorem 3 (ID consistency) Consider the PDE 9tU = Cu + g(u) for 

some local, isotropic, homogeneous, linear operator C, and for some smooth 
nonlinear reaction g . Model the dynamics on overlapping elements of an equi- 
spaced grid Xi = ih. Let Vi(x,t) denote the subgrid field in the ith element 
satisfying the PDE 9tVi = Cvi -\- g(vi) on the interval (Xt_i, Xi+i) with the 
moderated interelement coupling conditions 

Vi(Xt±i,t) =yvt±i(X|±i,t) + (l -y)vt(X|,t). (24) 

When interelement interactions are truncated to residuals O(y^) the grid 
values Ut(t) = Vi{Xi,t), at full coupling y — ] , evolve consistently with the 
PDE 9tU — Cu + g{u) . 

Proof: We proceed with some classic operator algebra pOl e.g.]. The 
principle obstacle is to transform subgrid spatial differences, indicated by 
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subscript x, into macroscale grid differences, indicated by subscript i. Begin 
with the PDE on the ith element: 9tVt = £vi + g(v|) . Because the operator £ 
is isotropic and homogeneous it may be formally expanded in even centred 
differences as 

oo 

k=0 

for some coefficients and corresponding function i. For example, in ap- 
plication to reaction diffusion PDEs, we would write the diffusion operator 
9^ = sinh^^ (-^Sx)] • Then provided the leading coefficient £2 7^ , the PDE 
9tVt = [£o + ^(Sx)Ni+ g("^i) is equivalently written [(9t — ^o)vi — g(vi.)] = 
S^Vi where l^^ is the inverse of function I (implicitly assumed invertible). 
Evaluate this last form of the PDE at x = so that the on the left- 
hand side becomes simply the grid value Ui and by the coupling condi- 
tions (Eljll the centred spatial difference on the right-hand side becomes 
the centred grid difference y6?Ui . This evaluation then gives the evolution 
l^^ [(9t — lo]W-i — g(U-i)] = yb\\li on the macroscale grid. 

Now reverting the inverse function, this grid evolution is equivalent to 



(25) 



k=0 



For example, for the diffusion equation 



9tU, = — [2sinh-^(^Vy5,)] Ut 



H2 
1 

h2 



2 3 4 



Thus a truncation of ( 1251) to errors C(yP) resuhs in a discrete model with 
stencil width of 2p — 1 . But specifically relevant to the theorem is that the 
equivalent differential equation of this discrete model evaluated at full cou- 
pling: the error in approximating C by the truncated version of (1231) (at full 
coupling y = 1) is dominated by the leading neglected term, namely lipb^ . 
As the element size h. — > 0, this error is 0(£2pH2^). For example, for the dif- 
fusion operator 9^, the coefficients ijk. = C(l /H^) and so the discrete model 
is consistent with the diffusion PDE to error O (h^^^^) as grid size h. — > . 



^If the leading coefRcient in the expansion of C is Ixn. 7^ , because the lower order 
coefficients are zero (or asymptotically small as in the Kuramoto-Sivashinsky pde), then 
more coupling conditions like (j24p couple with the next nearer neighbouring elements. 
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For interest, immediately generalise this theorem to a class of nonlinear 
reaction-dissipation pdes. 

Corollary 4 Consider the pde 9tU = £F(u)+g(u) for some local, isotropic, 
homogeneous, linear operator C, and for some smooth nonlinear F and reac- 
tion g . Model the dynamics on overlapping elements of an equi-spaced grid 
X| = ih. Lei V|(x, t) denote the subgrid field in the ith element satisfying the 
PDE 9tVt = £F(v|) + g[vi) on [Xi_-\,Xi^-\) with the moderated interelement 
coupling conditions 

F[vt(Xt±i,t)] =yF[v,±i(X,±i,t)] +(1 -y)F[vi(X,,t)] . (26) 

When interelement interactions are truncated to residuals O(y^) the grid 
values Ui(t) = vt(Xi,t), at full coupling ■y — ^ , evolve consistently with the 
PDE 9tU = CY[u) + g(u) . 

Proof: Naturally generalise the proof for Theorem [HI l|k 

In this theorem and corollary the subgrid microscale operator C need not 
be differential. For example, C could be a microscopic discretisation as in 
the numerical construction of the previous Section HI for example, 

^subgrid =4sinh^ [^HsubgridQx] ==4sinh^ 

The proof still holds. Further, the proof still holds if the subgrid field is not 
a scalar. This last observation empowers the following corollary on general 
consistency in two dimensions that Section [3] observed specifically for the 
Ginzburg-Landau PDE.. 

Corollary 5 (2D consistency) Consider a reaction- diffusion pde 9tU = 
V^u + g (u) ( such as the Cinzburg-Landau equation ([2]) ) modelled on over- 
lapping elements with subgrid fields v-tj(x,'y,t] coupled by conditions (jl])- 
1^. When the interactions are truncated to residuals O(y^) the grid values 
Utj(t) = vtj(Xi, Yj, t), at full coupling y = 1 , evolve consistently with the 
PDE to error O (jn^^^-^) . 

Proof: Apply the previous Theorem [3] twice. First, treat coordinate y 
as a parameter so that Iq — 9^ and ^[^^] — 9^. Then by Theorem [3] the 
semi-discrete 'grid' values Vij(Xt,-y,t), for discrete i and parametrised by 
continuous y, evolve consistently with the reaction-diffusion PDE. Second, 
treat index i as a parameter and consider the discrete modelling in coordi- 
nate y so that now £o involves operators acting on the x-grid and ^(6^) — ^u- 



2^' 



= 4sinh^ [^sinh^^ ifii] 
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Then by Theorem [3] the 2D grid values Utj = Vij(Xi,Yj,t) evolve consis- 
tently with the semi-discrete system generated in the first step, which in 
turn is consistent with the diffusion pde. The observation at the end of the 
proof for Theorem [3] provides the order of error. d|k 



6 Conclusion 

We explored novel macroscale discretisation of reaction-diffusion dynamics 
in two spatial dimensions. This work generalises considerable earlier work 
on modelling dynamics in one dimension. By dividing space into overlap- 
ping elements, we have shown that the specific interelement coupling condi- 
tions (jlj)-(l5l) have important theoretical and practical consequences. 

Section [2] discussed how these coupling conditions ensure that centre man- 
ifold theory applies to assure us of the existence, relevance and approximation 
of the slow manifold that is the macroscale discretisation of general reaction- 
diffusion equations in two dimensions. Further, Section proved that the 
resulting discrete models will also be consistent, as the macroscale grid size 
h — > , with the continuum or microscale dynamics. Thus the holistic dis- 
cretisations generated in this novel approach have the dual justification of 
both consistency and centre manifold theory. 

This strong theoretical support appears to be straightforwardly generalis- 
able to reaction-diffusion dynamics in three or more dimensions. The support 
also appears to be straightforwardly generalisable to higher order pdes, just 
as the theory supports the discrete modelling of one dimensional higher order 
PDEs such as the Kuramoto-Sivashinsky PDE [TU [16] . An interesting issue 
for further research is whether there are alternative interelement coupling 
conditions in two or more dimensions that have additional desirable prop- 
erties. For example, can one find interelement coupling that generates the 
so-called compact discrete approximations to the diffusion Laplace operator? 
rather than the spatially extended ones generated here. 

Sections [H] and HI using the example of the real Ginzburg-Landau pde, 
explored many of the technical issues necessary to apply the approach. In 
contrast to one spatial dimension, we found that a purely algebraic approach 
can only model to low order of accuracy. Although such low order models 
do reasonably accurately predict the dynamics, as seen in bifurcation dia- 
grams, higher order accuracy is desirable. Consequently we introduced and 
explored an approach where the microscale subgrid dynamics are described 
numerically, but with algebraic expressions for coefficients so that we con- 
struct an algebraic model for the macroscale discretisation. The computer 



27 



algebra package used does make a difference: we found reduce more than 
an order of magnitude faster than Mathematica. Further research may 
see if iterative techniques based upon the sparse microscale interactions are 
more effective than the direct LU factorisation employed here. This work 
provides a novel approach and theory for the sound and accurate closure of 
macroscale discretisations. 



Acknowledgement The Australian Research Council Discovery Project 
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